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Abstract 

Background: The emergence of pyrethroid resistance in the nnalaria vector, Anopheles arabiensis, threatens to 
undernnine the considerable gains nnade towards eliminating malaria on Zanzibar. Previously, resistance was 
restricted to the island of Pemba while mosquitoes from Unguja, the larger of the two islands of Zanzibar, were 
susceptible. Here, we characterised the mechanism(s) responsible for resistance on Zanzibar using a combination of 
gene expression and target-site mutation assays. 

Methods: WHO resistance bioassays were conducted using l-5d old adult Anoplieles gambiae s.l. collected between 
201 1 and 2013 across the archipelago. Synergist assays with the P450 inhibitor piperonyl-butoxide were performed 
in 2013. Members of the An. gambiae complex were PCR-identified and screened for target-site mutations 
(/cdr and Ace-1). Gene expression in pyrethroid resistant An. arabiensis from Pemba was analysed using 
whole-genome microarrays. 

Results: Pyrethroid resistance is now present across the entire Zanzibar archipelago. Survival to the pyrethroid 
lambda-cyhalothrin in bioassays conducted in 2013 was 23.5-54.3% on Unguja and 32.9-81.7% on Pemba. We 
present evidence that resistance is mediated, in part at least, by elevated P450 monoxygenases. Whole-genome 
microarray scans showed that the most enriched gene terms in resistant /\n. arabiensis from Pemba were associated 
with P450 activity and synergist assays with PBO completely restored susceptibility to pyrethroids in both islands. 
CYP4G16 was the most consistently over-expressed gene in resistant mosquitoes compared with two susceptible strains 
from Unguja and Dar es Salaam. Expression of this P450 is enriched in the abdomen and it is thought to play a role in 
hydrocarbon synthesis. Microarray and qPCR detected several additional genes putatively involved in this pathway 
enriched in the Pemba pyrethroid resistant population and we hypothesise that resistance may be, in part, related to 
alterations in the structure of the mosquito cuticle. None of the /cdr target-site mutations, associated with pyrethroid/DDT 
resistance in An. gambiae elsewhere in Africa, were found on the islands. 

Conclusion: The consequences of this resistance phenotype are discussed in relation to future vector control strategies 
on Zanzibar to support the ongoing malaria elimination efforts on the islands. 

Keyw/ords: Anoplieles arabiensis, Zanzibar, Pyrethroid resistance, P450s, Gene expression 



^ Correspondence: hranson@liverpool.ac.uk 

^Department of Vector Biology, Uverpool School of Tropical Medicine, 
Pembroke Place, Uverpool L3 5QA, UK 

Full list of author information is available at the end of the article 

O© 2013 Jones et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative 
BIoIVIGCI CGntrsI commons Attribution License (http://creativecommons.0rg/licenses/by/2.O), which permits unrestricted use, distribution, and 
reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication 
waiver (http://creativecommons.0rg/publicdomain/zero/l.O/) applies to the data made available in this article, unless otherwise 
stated. 



Jones et al. Parasites & Vectors 2013, 6:343 
http://www.parasitesandvectors.conn/content/6/1/343 



Page 2 of 1 3 



Background 

A major obstacle facing malaria control is pyrethroid re- 
sistance in the malaria vectors, Anopheles gambiae and 
Anopheles funestus [1,2]. Populations of highly resistant 
An, gambiae span West Africa [3,4] while pockets of resist- 
ance continue to emerge throughout the east of the con- 
tinent [5-7]. With alternative compounds for long-lasting 
insecticide nets (LLINs) and indoor residual spraying (IRS) 
still in the development pipeline, selection for resistance 
looks certain to continue unabated unless effective resist- 
ance management strategies are implemented. 

The search for genetic markers underpinning insecticide 
resistant phenotypes in natural insect populations poses 
significant challenges. This is exemplified in An. gambiae 
which displays a profile of resistant phenotypes across 
Sub-Saharan Africa. Identifying and tracking resistance- 
associated genetic markers is identified as a priority in the 
World Health Organisations (WHO) Global Plan for In- 
secticide Resistance Management (GPIRM) [1]. Resistance 
to pyrethroids and DDT is strongly associated with the 
knockdown resistance {kdr) target-site mutations (L1014F/ 
S) in the voltage gated sodium channel [8-10] but these 
mutations cannot account for all the variation observed in 
resistance bioassays and indeed, are absent from some re- 
sistant populations [11,12]. Pyrethroid resistance is medi- 
ated in part at least by the elevated expression of metabolic 
gene families (e.g. P450 monooxygenases/glutathione-S- 
transferases/esterases) and transcriptional profiling of re- 
sistant and susceptible An, gambiae followed by functional 
validation has begun to pinpoint individual genes respon- 
sible [13,14]. 

Substantial reductions in malaria have been made in 
Zanzibar through (i) the administration of artemisin com- 
bination therapies (ACTs), (ii) island wide distribution of 
LLINs and (iii) IRS campaigns using the pyrethroid, 
lambda-cyhalothrin [15,16]. In 2010, pyrethroid resistance 
was reported in Anopheles arabiensis from the Zanzibar 
archipelago [17] but curiously, resistance was confined to 
the smaller island of Pemba, which lies 50 km north of the 
larger island Unguja. It was uncertain whether the dis- 
crepancy in resistance levels between Pemba and Unguja 
was due to different selective forces acting on the islands. 

In this study, we set out to characterise the mechan- 
ism(s) responsible for the ongoing selection of resistance 
on Zanzibar using DNA and RNA-based approaches. 
Furthermore, we continued to monitor for phenotypic 
resistance across the islands between 2011 and 2013 in 
view of ongoing strategies for managing the emergence 
of pyrethroid resistance in Zanzibar. 

Methods 

Insect collections 

Anopheles larvae were collected from breeding pools 
across Unguja and Pemba during the main rainy season 



(April/May) between 2011 and 2013 (Figure 1). The sites 
were chosen to represent the typical diversity of breed- 
ing grounds present on the islands and their inclusion 
varied year to year depending on the number of larvae 
available. Four sites were sampled from Unguja: Mwera, 
Chuini, Kilimani, and Cheju. Breeding pools from 
Mwera and Chuini were generally from irrigated rice 
fields while Cheju is situated within the Jozani forest 
consisting of varied tropical vegetation. Kilimani is situ- 
ated on the coastal plain of the major capital Stone 
Town where the tide frequently submerges breeding 
sites resulting in brackish water. The breeding sites from 
Pemba were located in Pujini, Kiungoni, Uwandani, 
Mangwena and Tumbe. The majority of mosquitoes 
from these were sampled within rural inland settings 
with the exception of Tumbe which lies close to the 
north coast of the island. Mosquitoes were sampled from 
different pools to reduce sampling bias from families. 
Larvae were transported to insectaries at ZMCP and 
reared to adults under 27°C ± 2°C. 



Phenotypic bioassays 

Resistance assays were performed on one to five day old 
non blood-fed female adult An, gambiae using WHO 
susceptibility tests [18]. In 2011 and 2013, An, gambiae 
were exposed to a selection of the following insecticide- 
impregnated papers: lambda-cyhalothrin (0.05%), per- 
methrin (0.75%), DDT (4%) and bendiocarb (0.1%). In- 
sufficient numbers of larvae were collected from some 
sites to test all insecticides and in these sites we priori- 
tised lambda-cyhalothrin as this insecticide has been 
sprayed widely across Zanzibar as part of the IRS cam- 
paign since 2006 (ZMCP unpublished report, 2011). 
Control assays using non-treated papers were run 
alongside all tests. Mortality was scored 24 h post- 
exposure. In 2013, we performed additional synergism 
assays with the P450 inhibitor piperonyl butoxide 
(PBO). Mosquitoes were pre-exposed to PBO (4%) for 
one-hour and then transferred to the insecticide expos- 
ure tube for an additional hour. Mortality was scored 
24 h after exposure to insecticide. Control assays expos- 
ing mosquitoes to PBO only were run in each synergist 
assay. 

In 2012, a more comprehensive assessment of lambda- 
cyhalothrin resistance was performed to quantify the dif- 
ference in resistance between An, gambiae from Pemba 
and Unguja. The lethal concentration for 50% mortality 
(LC50) was determined from two sites on Unguja 
(Chuini & Mwera) and Pemba (Mangwena & Tumbe). 
Mosquitoes were exposed to insecticide impregnated 
papers at a range of five concentrations of lambda- 
cyhalothrin (0.005% - 0.2%) for one hour as described 
above. 
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An, gambiae complex identification and 
resistance-associated SNPs 

Members of the An, gambiae species complex were deter- 
mined using allele specific-PCR [19]. Target-site mutations 
in the sodium channel (resistance to pyrethroids/DDT) 
and acetylcholinesterase {ace-l) (resistance to carbamates/ 
organophosphates) were screened using TaqMan allelic 
discrimination assays [20,21]. 

Whole genome microarrays 

To investigate whether differential gene expression under- 
scored the resistant phenotype on Pemba, two whole gen- 
ome microarray experiments were performed in 2011 and 
2012. The emergence of pyrethroid resistance on Pemba, 
but apparent absence from Unguja, permitted a direct 
comparison of RNA extracted from mosquitoes from each 
island. However, in 2011 we were constrained by the num- 
ber of mosquitoes collected from Unguja and presence of 
An. gambiae s.s. from our main site Kilimani. Therefore, 
surviving mosquitoes exposed to lambda-cyhalothrin were 
compared directly with the insecticide susceptible An, 



arabiensis lab-colony from Mozambique (MOZ) [11]. In 
2012, the experimental design was extended to include 
mosquitoes from (i) Pemba (Mangwena), (iii) Unguja 
(Chuini) and (iii) a susceptible field strain of An, arabien- 
sis from Dar es Salaam (DAR). Four day old, non blood- 
fed An, arabiensis from Ilala and Kinondoni in Dar es 
Salaam were stored in RNAlater® 24 h post-exposure to 
insecticide and display no resistance to pyrethroids or 
DDT. The inclusion of DAR provided an independent sus- 
ceptible population to improve the strength of resistant 
gene expression association. Three biological replicates 
from each strain were integrated into an interwoven loop 
design as this has been shown to provide robust power in 
detecting differences in gene expression [22]. 

The An, gambiae 15 k whole genome microarray was 
used in all experiments [14]. Total RNA was extracted 
from batches of 8-10 female mosquitoes using the 
Ambion RNA4PCR kit. RNA was labelled with cy3 and 
cy5 dyes and hybridised to the microarray chip accord- 
ing to the protocol described previously [14]. RNA and 
cRNA quantity and quality were assessed using a Nano- 
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drop spectrophotometer (Nano-drop Technologies) and 
Bioanalyser (Agilent) respectively. 

Microarray analysis 

The direct comparison between PEM and MOZ was 
analysed using GeneSpring v. 11.5 (Agilent) which applies 
a standard t-test to normalised raw fluorescence values. 
Following normalisation with LIMMA [23] the MAA- 
NOVA package in R [24] was applied to the loop design 
as used previously [14]. The two independent experi- 
ments from 2011 and 2012 were considered together, 
providing independent datasets from which improved 
confidence in identifying genuine candidate genes could 
be made. Probes were filtered based on the hypothesis of 
significance (q<0.05) and under the hypothesis that ex- 
pression should be consistently greater in the resistant 
compared to the susceptible populations (i.e. PEM > UNG/ 
DAR). Microarray data are available in the ArrayExpress 
database (www.ebi.ac.uk/arrayexpress) under accession 
number E-MTAB-2075 (PEM vs. MOZ direct comparison) 
and E-MTAB-2074 (PEM vs. UNG/DAR loop design). 

Patterns of resistance-associated gene expression on 
Zanzibar 

The expression levels of candidate genes for pyrethroid 
resistance from the microarray analysis were assessed 
using quantitative reverse transcription PGR (qRT-PCR) 
following the MIQE guidelines [25]. 

qPCR was performed on mosquitoes collected from 
2012 and 2013 with three objectives: 

I. To independently validate gene expression in An. 
arabiensis from 2012 in the microarray experiments. 

II. To determine gene expression patterns of candidate 
genes from Pemba and Unguja in collections made 
in 2013. 

IILTo independently validate the expression of a subset 
of genes up-regulated in Pemba which are putatively 
part of the hydrocarbon synthesis pathway. 

Whole female mosquitoes were collected 24 h post- 
exposure to insecticide or non-treated papers, briefly 
chilled (10-20 min) and stored in RNAlater® according 
to the manufacturers instructions. Total RNA was ex- 
tracted from An, arabiensis using PicoPure® RNA Isola- 
tion Kit (Invitrogen) and treated with DNase I (Qiagen). 
Three biological replicates of RNA extracted from ten 
mosquitoes were used for the 2012 samples to replicate 
the microarray experiments. The following year, four- 
five biological replicates of five mosquitoes per RNA 
sample were used. RNA quantity and quality was 
checked using a Nano-Drop spectrophotometer (Nano- 
drop Technologies) and Bioanalyser (Agilent) respect- 
ively. cDNA was synthesised from -0.5-1 (ig of RNA 



using oligo(dT)2o (50 \M) and Superscript III (200 U) 
(Invitrogen) and purified through a DNA-binding col- 
umn (Qiagen). 

Multiple pairs of exon-exon spanning primers for tar- 
get and control genes were designed in silico against the 
Anopheles gambiae PEST sequence (Taxon ID 180454) 
using primer-BLAST (NCBI). The PGR efficiency, dy- 
namic range and specificity of the primer pair were cal- 
culated from running a standard curve over a five-fold 
dilution series (input cDNA from x5'^ to x5'^). The pri- 
mer sequences, amplicon length, efficiencies and gene 
location are given in Additional file 1. 

PGR reactions (20 (il) were performed on the MXPro 
qPGR system (Agilent) with 10 \A Brilliant III SYBR 
Green (Agilent), 300 nM of primers, 2.5 (il of input 
cDNA (diluted 100-fold) and the total volume made up 
with sterile-distilled water. The thermal profile used 
throughout consisted of 95°G for 10 min followed by 
40 cycles of 95°G for 10s and 60°G for 10s. Melt curves 
were run after each PGR to ensure the specificity of the 
amplified products. Three technical replicates were run 
for each sample and no template controls were run for 
each gene. 

Data analysis of gene expression 

The data were pre-processed prior to analysis. A Gq 
value of 35 was considered our limit of detection (LOD) 
and any samples that failed to amplify above this thresh- 
old were given a value of 35. Outliers were defined as 
those values which give a standard deviation above 0.5 
from the three technical replicates and were removed 
from the dataset. 

Once pre-processed the Gq values were adjusted accord- 
ing to the efficiency of the primer pair. The qPGR repeats 
were averaged and data normalised against the average 
values for the ribosomal S7 protein (AGAP010592) and 
ubiquitin (AGAP007927). The relative quantities of each 
sample were calculated and log-transformed to normalise 
the distribution of the data for parametric statistical ana- 
lysis [26]. A one-sided t-test at p < 0.05 was used to assign 
significance between treatments. 

Copy number analysis of P450 candidates 

Several insects have adapted to selection from insecti- 
cides by evolving variations in gene copy number (other- 
wise known as copy number variation (GNV)) (reviewed 
in [27]). Quantitative PGR (qPGR) using SYBR Green 
was used to determine whether the elevated expression 
of CYP4G16, CYP6Z2, and CYP6Z3 in An, arabiensis 
from Pemba is due to gene amplification. The quantity 
of genomic DNA (gDNA) extracted from individual An, 
arabiensis from Unguja and Pemba was analysed using a 
PicoGreen® assay [28] and diluted to 1 ng/(il. Internal 
exonic primers were designed for the target genes 
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CYP4G16, CYP6Z2 and CYP6Z3 and for two inter- 
nal reference genes, elongation factor (AGAP005128) 
and glucose-6-phosphate dehydrogenase (AGAP010739) 
(Additional file 1). The efficiency, dynamic range and spe- 
cificity of each primer pair were assessed by performing a 
standard curve on a dilution series of DNA (five-fold dilu- 
tions from x5'^ to x5'^). SYBR Green qPCR reactions were 
performed on the Stratagene MXPro using the exact same 
conditions as RT-qPCR (see above) with the exception 
that 2.5 [A input gDNA was added to each reaction. To 
improve the power of detecting a real difference in copy 
number between resistant and susceptible mosquitoes, 
copy number analysis was performed on eight individual 
An. arabiensis which died at the lowest concentrations of 
insecticide exposure (0.001-0.005%% \-cyhalothrin) from 
Unguja and ten individuals which survived the highest 
concentrations (0.1-0.2% \-cyhalothrin) from Pemba. The 
average gene copy number was calculated using the AACq 
method described previously [29]. 

Sodium channel sequencing 

Exon 20 of the sodium channel harbours key residues 
for pyrethroid and DDT binding in the so-called 'binding 
pocket' and mutations in this region confer resistance in 
other insect pests [30]. To determine whether pyrethroid 
resistance was due to previously unidentified SNPs in 
this region of An, arabiensis^ a -400 bp fragment of 
Exon20 was amplified from gDNA extracted from indi- 
viduals on Pemba and Unguja collected in 2012. PGR 
reactions were undertaken in 25 \A total volumes 



consisting of 0.5U Taq Polymerase (KappaBiosystems), 
0.5 \A dNTPs (10 mM) and 400 nM of forward and re- 
verse primers. Two alternative splice forms of this exon 
occur in An, gamhiae (Exon c/d) and two sets of primer 
pairs were designed in flanking introns to amplify each 
variant (Additional file 1). Thermal cycling parameters 
were 95°C for 2 min followed by 40 cycles of 95°C for 
30 s, 55°C for 30 s, and 72°C for 30 s and a final extension 
step of 72°C for 5 mins. Samples for sequencing were 
performed by Macrogen (Amsterdam) with the forward 
and reverse primers used as the sequencing primers. 

CYP4G16 cloning and sequencing 

The full length CYP4G16 was amplified from cDNA syn- 
thesised from pools of five An, arabiensis individuals from 
Pujini (Pemba) and Mwera (Unguja). The PGR products 
were amplified using the forward and reverse primers 
AGAP001076_F1/R1 (Additional file 1) and cloned into 
the pJET 1.2/blunt cloning vector (Fermentas). Two 
clones from each strain were sequenced (Gore Genomic 
Facility, University of Sheffield) with the pJET forward and 
reverse primers as well as three nested forward primers 
(AGAP001076_F2/F3/F4) (Additional file 1). 

Results 

Study area and current distribution of >An. gambiae 
complex on Zanzibar 

According to the sites we visited during the study 
(Figure 1) An, arabiensis is the dominant vector on Zanzibar 
(Table 1). The exception was in Kilimani, situated within 



Table 1 Species identification of a subset of ^n. gambiae s.l. collected between 2011 and 2013 on Zanzibar 



Year 


Island 


Site 


No. 


An. 


%An. 


An. 


%An. 


An. 


%An. 


An. 


%An. 








tested 


arabiensis 


arabiensis 


merus 


merus 


gambiae 


gambiae 


quad^ 


quad^ 


2011 


Unguja 


Kilimani 


78 


13 


16.7 


12 


15.4 


53 


67.9 


0 


0.0 




Pemba 


Cliwale 


24 


23 


95.8 


0 


0.0 


0 


0.0 


0 


0.0 






Kiungoni 


24 


24 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 






Uwandani 


24 


24 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 


2012 


Unguja 


Cliuini 


72 


72 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 






Mwera 


110 


110 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 




Pemba 


Tumbe 


340 


206 


60.6 


134 


39.4 


0 


0.0 


0 


0.0 






Mangwena 


286 


284 


99.3 


4 


1.4 


0 


0.0 


0 


0.0 


2013 


Unguja 


Cliuini 


102 


102 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 






Mwera 


71 


71 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 






Kilimani 


66 


66 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 






Cheju 


28 


16 


57.1 


0 


0.0 


6 


21.4 


6 


21.4 




Pemba 


Pujini 


75 


75 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 






Kiungoni 


92 


92 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 






Uwandani 


69 


69 


100.0 


0 


0.0 


0 


0.0 


0 


0.0 






Tumbe 


71 


70 


98.6 


0 


0.0 


1 


1.4 


0 


0.0 



^Anopheles quadriannulatus. 
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the main urban area of Stone Town and adjacent to the 
coast, where we identified 68% An. gambiae s,s,, 17% An, 
arabiensis and 15% An, merus in 2011. Two years later all 
mosquitoes were scored as An, arabiensis from this site 
(N = 66). In Tumbe, on the north coast of Pemba, 39% of 
Anopheles were typed as An, merus (N = 134) consistent 
with their salt-water adaptation. However, the following 
year, we did not identify a single An, merus from Tumbe 
(N = 71); 98% samples from this site were An, arabiensis, 
2% An, gambiae. 



Insecticide resistance and synergism 

Following the emergence of pyrethroid resistance on 
Pemba [17] we generated lambda-cyhalothrin LC50 
curves from two sites on each island to quantify the 
scale of the resistance phenotype between the islands 
(Figure 2; Additional file 2). The LC50S from Pemba 
((Mangwena = 0.121% (95% CI = 0.080-0.184) and Tumbe = 
0.092% (95% CI = 0.060-0.141%)) were a magnitude 
greater than those from Unguja ((Mwera = 0.007% (95% 
CI = 0.005-0.010%) and Chuini = 0.009% (95% CI = 0.007- 




Mangwena, Pemba Tumbe, Pemba 




0.001 0.005 0.020 0.050 0.200 0.001 0.005 0.020 0.050 0.200 



Log Concentration (%) Log Concentration {%) 

Figure 2 Dose-response curves for An. arabiensis exposed to lambda-cyhalothrin on Zanzibar. An. gombioe were exposed to a range of 
concentrations of lambda-cyhalothrin in WHO susceptibility tests. The lethal concentration for 50% mortality (LC50) was calculated by fitting a 
generalised linear model (GLM) using a binomial logit-link function. The LC50 is shown above each plot. For Tumbe, two additional curves are 
shown specific for An. arabiensis (red dashed) and An. merus (orange dashed). 

V ) 
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0.011%)) giving resistance ratios of between 10 and 17-fold 
A re-assessment of the data from Tumbe showed that An. 
merus were markedly more sensitive to lambda-cyhalothrin 
displaying an LC50 closer to that of An. arabiensis from 
Unguja (LC50 = 0.010%; 95% CI = 0.006-0.018%) (Figure 2). 

In addition to the dose-response curves 1125 An. 
gambiae were tested against the WHO diagnostic dose 
of lambda-cyhalothrin (0.05%) in susceptibility bioassays 
in 2011 and 2013 (Figure 3; Additional file 2). WHO- 
defined resistance (less than 90% mortality [18]) was ob- 
served throughout the study in all sites from Pemba 
(Figure 3). The lowest mortalities observed were from 
Chwale in 2011 (19.0%; 95% CI = 12.1-28.3) and Pujini 
in 2013 (18.3%; 95% CI 11.3-27.9%). Although these data 
indicate that resistance remains particularly strong on 
Pemba, significantly higher mortality levels were observed 
in three out of the four sites tested on Pemba in 2013 
compared with Pujini (x^ p < 0.001 for each pair- wise 
comparison). This suggests that resistance may not be 
homogeneously expressed across the island. In Kiungoni, 
mortality to lambda-cyhalothrin approximately doubled 
between 2011 and 2013 indicating that in this site at least 
that resistance decreased over the three year study period. 

As reported earlier [17], An. arabiensis from Unguja 
were susceptible to pyrethroids at the start of this study 
(2011). However, resistance had emerged on the island 
(81.9% - 88.4% mortality) in 2012. In the subsequent 
year, further reductions in mortality were observed in 
Chuini (45.7%; 95% CI 36.1-55.7%) and Mwera (67.0%; 
95% CI 57.1-75.6%) confirming that resistance had arisen 
during the course of the study. 

Synergist assays with piperonyl butoxide (PBO 4%), a 
general inhibitor of P450 monooxygenases, completely 
restored susceptibility to lambda-cyhalothrin in Unguja 
and Pemba (mortality = 96-100%) (Figure 3). This sug- 
gests that P450s may contribute towards the resistance 
phenotype we observe in An. arabiensis. Unfortunately, 



due to a restricted number of insects available, we were 
only able to perform one set of assays using DDT but 
the 2011 data do suggest potential cross-resistance in 
Pemba (Kiungoni, 78.8% mortality). 

Finally, the carbamate bendiocarb has recently re- 
placed lambda-cyhalothrin in IRS on Zanzibar [17] and 
our bioassay data show that all mosquitoes tested to date 
are completely susceptible to this insecticide (100% mor- 
tality, Unguja N = 57; Pemba N = 192; Additional file 2). 

Existing target-site markers for pyrethroid resistance 

A sub-sample of Anopheles were screened for the kdr 
target site mutations (1014 F and 1014S) and ace-lR. 
None of these known resistance alleles were found over 
the three year period (Unguja N = 120; Pemba N = 251). 
Sequencing exon 20 of the sodium channel which harbours 
target-site mutations in other insects [30] failed to identify 
any additional non-synonymous mutations suggesting that 
alternative mechanisms underlie the pyrethroid resistance 
phenotype in An. arabiensis from Zanzibar. 

Gene expression in An. arabiensis from Pemba Island 

A direct comparison of the adult transcriptome between 
An. arabiensis collected from Pemba in 2011 and a sus- 
ceptible laboratory colony originating from Mozambique 
(MOZ) yielded 2214 probes significantly differentially 
expressed (1071 up (48.3%); 1143 down (51.7%)) 
(Additional file 3). This large number of probes and 
the extremely high fold changes at the ends of the 
expression distribution most probably reflect genetic 
divergence between the field and laboratory strains. 
We therefore do not consider this experiment in isolation; 
however, under the assumption that the mechanism re- 
sponsible for pyrethroid resistance was the same in 2011 
and 2012, this data set was combined with a second 
experiment conducted in 2012 on An. arabiensis from 
Pemba (PEM), Unguja (UNG) and pyrethroid-susceptible 
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Figure 3 Susceptibility of An. gambiae to lambda-cyhalothrin in Zanzibar between 201 1 and 2013. Mosquitoes were exposed to 
lambda-cyhalothrin (0.05%) in WHO susceptibility tests between 201 1 and 2013. Blue (Unguja) and red (Pemba) percentage mortalities with 
95% binomial confidence intervals are given one-hour of exposure to insecticide. Darker bars represent synergist assays using one-hour 
pre-exposure with piperonyl butoxide (PBO) performed in 2013. 
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mosquitoes from Dar es Salaam (DAR). Using the 2214 
probes from the 2011 experiment as our initial probe list, 
we subsequently filtered probes based on the hypothesis 
of significance in each pair-wise comparison (false discov- 
ery rate (FDR) adjusted p-value < 0.05) and greater expres- 
sion in the resistant strain PEM compared to UNG and 
DAR (N = 2645). This left just 268 probes representing 
208 unique transcripts (Additional file 3). 

The most over-expressed transcript in the PEM vs 
DAR comparison was CYP4G16 (AGAP001076-RA) 
(fold-change (PC) = 5.4). CYP4G16 has four alternatively 
spliced transcripts and probes specific for each of these 
were among those most highly expressed in PEM against 
DAR (average PC = 4.4) and UNG (PC = 2.0) (Pigure 4). 

The most over-expressed transcripts in PEM compared to 
UNG included a Niemann-pick type c gene (AGAP002852- 
RA), two transcripts for enzymes involved with 5' nucleotid- 
ase activity (AGAP005457-RA & AGAP005458-RA) and 
the alkaline phosphatase (AGAP001684-RA) and the P450 
CYP6Z2 (AGAP008218-RA) (Additional file 3). Two mem- 
bers of the CYP6Z family of P450s were overexpressed in 
all three PEM comparison and showed higher expression 
against UNG than samples from DAR; CYP6Z2 (PC = 2.6 
vs. 1.6) and CYP6Z3 (PC = 1.9 vs. 1.6) (AGAP008217-RA). 
Members of this class of CYP genes have been associated 
with pyrethroid resistance previously [31,32]. 

Functional annotation of the 268 candidate probes 
with the DAVID analysis software tool [33] yielded no 
significantly enriched biological terms. Under the stringent 
filtering strategy used here there is a strong possibility of 
missing genes or biological processes involved with the 



resistance phenotype. Therefore, we ran a DAVID ana- 
lysis on common probes over-expressed in the resistant 
mosquitoes from the 2012 experiment only (N = 2645) 
(Additional file 3). Following multiple hypothesis testing 
the most enriched cluster of genes was associated with 
P450 metabolism (enrichment score (ES) = 4.55), mito- 
chondrial processes (ES = 3.93) and ribosomal genes 
(ES = 3.36) (Additional file 4). A look at individual genes 
solely from the 2012 analysis found that no additional 
P450s were contained within the PEM over expressed 
subset than those already described above. 

Validation of candidate gene expression in Zanzibar 

Independent validation of expression data from high- 
throughput transcriptomic studies is a necessary step to- 
wards identifying and verifying candidate genes. Based 
on our final list of candidates we took six target genes 
forward for reverse-transcription quantitative PCR (RT- 
qPCR); the P450s CYP4G16, CYP6Z2 and CYP6Z3 along 
with three genes with putative involvement as part of a 
hydrocarbon synthesis pathway, acyl-coa thioesterase 
(AGAP003848), acyl-coa dehydrogenase (AGAP005662) 
and 3-hydroxyacyl-coa dehydrogenase (AGAP007784). 

RT-qPCR on An. arabiensis collected in 2012 from 
Mangwena (the same site used in the microarray analysis) 
confirmed a significant over-expression of CYP4G16 
in Pemba compared to Unguja (ddCt = 5.5; p < 0.014) 
but could not validate the over-expression of CYP6Z2 
(ddCt = 1.5; p = 0.127) or CYP6Z3 (ddCt = 0.9; p =0.206) 
(Figure 5A). 
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Figure 4 Volcano plots for expression between An. arabiensis from Pemba against samples from Dar es Salaam and Unguja. Expression 
plots for the 268 candidate probes significantly up-regulated in (i) An. arabiensis from Pemba collected in 2012 against UNG and (ii) DAR. Probes 
for the three candidate P450s are highlighted on each plot. 
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Figure 5 Quantitative PGR validation of candidate P450s in An. arabiensis from Zanzibar. A) Quantitative PCR (qPCR) analysis of CYP4G16, 
CYP6Z2 and CYP6Z3 was performed on An. orobiensis collected from Pemba (Mangwena) and Unguja (Chuini & Mwera) in 2012. The mean ± SEM 
for three ddCq values relative to Unguja are presented. (B) CYP4G16 expression in An. arabiensis collected in 2013. Four groups of mosquitoes 
were included in the analysis: (i) Pujini exposed to lambda-cyhalothrin (0.05%) (PujiniA) (ii) Pujini unexposed to insecticide (Pujini_untreated) (iii) 
Mwera (iv) Chuini. The mean ddCq values ± SEM of five biological replicates are presented. NS = non-significant, ""p <0.05 one-sided 
Student's t-test. 



In 2013, we collected mosquitoes in RNAlater from 
two sites on Unguja (Mwera and Chuini) and from the 
most resistant site on Pemba, Pujini. When we compare 
expression against An. arabiensis from Mwera - the least 
resistant site in 2013 (mortality to lambda-cyhalothrin 
67.0%) - RT-qPCR estimated a 1.6- and 1.5-fold increase in 
CYP4G16 in lambda-cyhalothrin-treated and non-treated 
mosquitoes respectively although neither difference was 
significant (p > 0.05) (Figure 5B). This is approximately 
a quarter of the difference in expression observed be- 
tween samples taken from Pemba and Unguja in 2012 
(Figure 5A) and could represent the increase of pyreth- 
roid resistance in Unguja by our sampling period in 
2013. 

Three genes with putative involvement in hydrocarbon 
synthesis which passed our filtering criteria in the micro- 
array analysis (acyl-coa thioesterase, acyl-coa dehydrogen- 
ase, hydroxyacyl-coa dehydrogenase) were analysed by 
qPCR on An. arabiensis collected in 2013. The patterns of 
up-regulation in Pemba were extremely similar between 
those detected through the microarray in 2012 and those 
analysed by qPCR in 2013 (Figure 6). This provides fur- 
ther independent validation of the array experiments and 
suggests that these genes could be co-expressed with 
CYP4G16 as part of the hydrocarbon synthesis pathway. 

Copy number variation of P450 candidates 

Copy number analysis of the three major candidate 
genes CYP6Z2, CYP6Z3 and CYP4G16 was performed 
on An. arabiensis from Pemba surviving the higher 
doses of lambda-cyhalothrin compared to susceptible 
counterparts from Unguja. None of the candidate genes 
showed any significant difference between the two 
groups and so the increase in expression observed in the 



microarrays is unlikely explained by copy number 
changes (Additional file 5). 

Discussion 

While national malaria control programmes targeting in- 
door resting/biting mosquitoes have enjoyed success at re- 
ducing malaria in recent years, this success has brought 
with it many challenges for the future of vector control in- 
cluding: a) the selection of physiological resistance [2], b) 
a change in Anopheles biting and resting habits [34] and c) 
shifts in the composition of vector species [35]. The evi- 
dence presented previously [17] and in this study suggests 
that at least two of these are threatening malaria control 
in Zanzibar. 

Anopheles arabiensis has supplanted An. gambiae s.s. as 
the major vector in many parts of East Africa where the 
two species are sympatric [36,37]. In the 1960s, it was re- 
ported that An. arabiensis had replaced An. gambiae s.s. 
as the major vector on Zanzibar following DDT spraying 
as part of the Global Malaria Eradication Programme [38] 
although by the 1990s, An. gambiae s.s. had re-established 
itself as the main vector on Zanzibar [39]. The data pre- 
sented here and in a previous publication [17] suggest that 
history is repeating itself, with An. arabiensis remerging as 
the dominant vector following the implementation of a 
scaled-up vector control programme. Pockets of An. gam- 
biae s.s. and An. merus were found in sites proximate to 
the coast (e.g. Kilimani (Unguja) and Tumbe (Pemba)) but 
only a few An. gambiae (N = 7) and no An. merus were 
identified in 2013. It is worth noting that our sampling 
strategy was restricted to collecting larvae from breeding 
pools and it important to broaden these collections to 
collect adult vectors indoors and outdoors, and screen 
these for the presence of Plasmodium before drawing 
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samples. For the 2013 qPCR data, the fold-changes are calculated 
from ddCq(Pemba)/ddCq(Unguja). 



conclusions on the role of alternative vectors in contrib- 
uting to the residual transmission on Zanzibar. 

Pyrethroid resistance is now present throughout the en- 
tire Zanzibar archipelago. In a previous study conducted 
in 2010-2011, resistance was detected solely in An. ara- 
biensis from Pemba suggesting that differential selective 
forces were acting on the two islands. However, we have 
subsequently shown that the proportion of mosquitoes 
surviving exposure to lambda-cyhalothrin (0.05%) in- 
creased approximately 3-fold in Unguja between 2012 
and 2013. This rise within a single year highlights how 
swift the spread and selection of resistance can occur. 

At present, it is difficult to determine the origins of re- 
sistance on Zanzibar. Pemba and Unguja are separated 
by a 50 km strip of the Indian Ocean and large numbers 
of human traffic moving between the two islands pro- 
vide a passive route for migration of anopheline mosqui- 
toes. Similarly, it is feasible that resistance could have 
emerged independently on Unguja. Resolving the genetic 
differentiation of Anopheles between Unguja and Pemba 
would not only provide information on the spread of re- 
sistance traits but present estimates of the diversity and 
effective population sizes as measures of vector control 
success. Regardless of the origin, resistance management 
strategies should be implemented equally on both islands 
to curb further selection. 

In bioassays conducted on Pemba we observed WHO- 
defined resistance (less than 90% mortality) in all sites 
tested. Nevertheless, with the exception of Pujini it was 
noticeable that average resistance levels dropped 
between 2011 and 2013 with three out of four sites 



showing mortality over 67% and large increases in pyr- 
ethroid mortality in some sites (e.g. mortality in An. ara- 
biensis to 0.05% lambda-cyhalothrin was 5.2% in Tumbe 
in 2012 but increased to 77.1% the following year). The 
vast majority of insects tested in these sites were An. 
arabiensis and therefore any increase cannot be attrib- 
uted to species differences. It is also difficult to ascribe 
the differences in resistance levels across Pemba to ex- 
perimental variation as we were careful to record the 
age of mosquito, as well as room temperature and hu- 
midity throughout - factors which can greatly influence 
mortality in bioassays [18]. 

It is possible that the change in insecticide use in IRS, 
from pyrethroids to carbamates may have contributed to 
the reduction in pyrethroid resistance. However, further 
rounds of bioassays, and an analysis of the strength of 
the resistance, via repeat determination of the LC50, is 
needed before it can be concluded that pyrethroid resist- 
ance has declined since the introduction of bendiocarb 
for IRS. 

In response to the discovery of pyrethroid resistance 
on Zanzibar we set out to identify the underlying mech- 
anisms in An. arabiensis to support future resistance 
management strategies. We performed two separate 
microarray experiments using stringent criteria to filter 
out probes based on significant over-expression in An. 
arabiensis from Pemba. This yielded just 208 transcripts 
and of these, replicate probes for three specific P450s 
{CYP6Z2, CYP6Z3 and CYP4G16) previously implicated 
in other resistance studies in An. gambiae were among 
the most over-transcribed genes. 

The P450, CYP4G16 was up-regulated in resistant An. 
arabiensis approximately 4.5- and 2.0-fold compared to 
DAR and UNG respectively. Although these are fairly 
modest changes in expression they were among the high- 
est observed in the comparison of the field populations. 
qPCR analysis gave an estimated 6-fold difference between 
An. arabiensis from Pemba compared to Unguja. Elevated 
expression levels of this enzyme were observed in families 
of An. arabiensis displaying an increased tolerance to del- 
tamethrin in Northern Cameroon [40]. A modest fold- 
change increase was observed in a laboratory resistant 
strain of An. arabiensis from Sudan using qPCR although 
this was not backed-up by the microarray data in the same 
study [41]. 

CYP4G16 is an alternatively spliced gene with four 
transcripts located on the X chromosome. Sequencing 
cloned CYP4G16 from pooled cDNA from Pemba and 
Unguja revealed no SNPs between the populations and 
showed that CYP4G16 is an extremely conserved en- 
zyme between An. arabiensis and An. gambiae. The ma- 
jority of insects only possess two CYP4G genes [42]. The 
other CYP4G gene in An. gambiae, CYP4G17, was also 
up-regulated in resistant An. arabiensis from Pemba 
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(1.9-fold versus DAR; 1.5-fold versus UNG) but was not 
significantly expressed in the comparison against MOZ 
in 2011. CYP4G genes have played a key evolutionary 
role in the terrestrial adaptation of Insecta by catalysing 
the final step of cuticular hydrocarbon synthesis and 
thus providing a protective waterproof layer [42]. The 
orthologue of CYP4G17 in Drosophila melanogaster 
(DmCYP4Gl) catalyses an oxidative decarbonylation of 
long-chain aldehydes, which serves as the final step of 
hydrocarbon synthesis from long-chain fatty acids in in- 
sects [43]. Synthesis of cuticular hydrocarbons from long- 
chain fatty acids is complex and requires a suite of elon- 
gases, reductases and dehydrogenases [44]. It is therefore 
of interest that the final list of candidates from the micro- 
array analysis included several genes involved in fatty-acid 
metabolism; acyl-coa thioesterase, acyl-coa dehydrogen- 
ase, hydroxyacyl-coa dehyrogenase, and over expression of 
these was confirmed in An. arabiensis collected in 2013 by 
qPCR. DmCYP4Gl is highly expressed in the fat body and 
carcass of male and female D. melanogaster whereas the 
orthologue of CYP4G16, DmCYP4GlS, is expressed in the 
thoraic- abdominal ganglion [45]. Tissue specific expres- 
sion in An. gambiae shows that both CYP4G genes are 
over expressed in the insect abdomen consistent with 
hydrocarbon synthesis in oenocytes (Ranson et al. unpub- 
lished; MozAtlas [46]). In the German cockroach, Blatella 
germanica, CYP4G19 was over-expressed five-fold in a 
pyrethroid resistant strain and expression was found to be 
greater in the abdomen than the head or thorax [47]. 
Cuticular-based resistance has been reported from several 
agricultural and medically important insects [48-50] with 
the over-expression of specific genes implicated in some 
cases [51,52]. However, whilst the up-regulation of CYP4G16 
does suggest a potential role for cuticular based resistance 
in An. arabiensis in Zanzibar, it is important to recognise 
that the functional work to validate this hypothesis is still 
underway. Whether or not PBO can synergise pyrethroids 
(and other insecticides) by inhibiting P450s essential to 
cuticular formation also needs further investigation. 

In contrast to target site resistance, where a single mu- 
tation can confer high levels of resistance [10] there is 
increasing evidence that insecticide resistance in Anoph- 
eles may not be underpinned by one single gene but by a 
suite of co-regulated enzymes in a more intricate path- 
way depending on the population. One example is the 
CYP6Z family in which members are over-expressed in 
pyrethroid resistant An. gambiae [31,32]. Functional work 
in these enzymes has found that they are unable to me- 
tabolise pyrethroids directly but instead, that they play a 
role in the clearance of pyrethroid metabolites. We 
propose that the coordinated up-regulation of multiple 
genes involved in CHC biosynthesis is a putative resist- 
ance mechanism underpinning the pyrethroid resistance 
in An. arabiensis in Zanzibar. 



The association of P450s with resistance on Zanzibar 
is supported by the near full restoration of pyrethroid 
susceptibility we observed in the PBO assays. PBO has 
been incorporated into a new generation of LLINs 
alongside pyrethroids to combat resistance [53]. While a 
wider adoption of these nets as IRM tools is restricted at 
present, given the resistance phenotype described in An. 
arabiensis in this study (i.e. over-expression of P450s 
with no kdr), it would be interesting to determine their 
effectiveness against wild caught mosquitoes on Zanzibar 
as a possible option for resistance management. 

Pyrethroid resistance in An. arabiensis remains patchy 
in both its distribution and severity in East Africa. A 
common phenotype in An. arabiensis of resistance to 
pyrethroids with no cross-resistance to DDT has now 
appeared in southern Uganda on the shores of Lake 
Victoria [54], in the Northern highlands of lower Moshi 
in mainland Tanzania [7]. There is evidence for a similar 
but milder resistant phenotype from Western Kenya 
[55]. The presence of DDT resistance on Pemba in 2011 
suggests that the resistance phenotype on Zanzibar dif- 
fers to that on the mainland. An additional bioassay con- 
ducted during the dry season (October) on Pemba in 
2013 confirms this finding (63.0% mortality; Additional 
file 2) and indicates that DDT resistance has remained 
relatively stable between 2011 and 2013. At present 
there is no evidence to suggest that kdr plays a role in 
any of these regions and in this study, a large screen for 
1014 F and 1014S (N = 371) as well as sequencing regions 
of the sodium channel where alternative kdr variants exist 
in other insect pests, provided no evidence to suggest the 
contrary. Whether a single mechanism of resistance in 
An. arabiensis is sweeping across the East African region 
or independent selection events are taking place is cur- 
rently unknown. Future population genetic approaches 
in An. arabiensis from Zanzibar and a wider geograph- 
ical region will no doubt provide clues on the spread of 
this resistance phenotype. 

Conclusion 

Insecticide resistance in malaria vectors is seen as a sig- 
nificant hurdle to malaria elimination in Africa. There 
are a few well cited examples of operational failure as a 
consequence of resistance and experimental evidence is 
gathering to show that LLINs and/or IRS are not ad- 
equately killing resistant vectors [56,57]. Haji et al. 
(2013) and colleagues have already shown that resistance 
has a significant operational impact on the effectiveness 
of LLINs and IRS and warned about its potential nega- 
tive effects on Zanzibar highlighting the importance of 
resistance management. Important proactive steps have 
been swiftly taken to address this issue on Zanzibar with 
the switch to alternative modes of action for targeted 
IRS. The plasticity in feeding and resting behaviour of 



Jones et al. Parasites & Vectors 2013, 6:343 
http://www.parasitesandvectors.conn/content/6/1/343 



Page 12 of 13 



An. arabiensis [58] will necessitate additional tools to 
control this vector. The case of Zanzibar represents an 
ideal opportunity to systematically monitor resistance, 
using entomological data and candidate molecular markers 
(described in this paper), and quantify its impact on our 
ability to control mosquito vectors. 
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